#include <gmp.h>
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <pthread.h>

/* It should be enough :D */
#define DEFAULT_N_ITERATIONS 25
#define DEFAULT_PRECISION    35000000 /* in bits */

#define PARAM_A 0
#define PARAM_B 1
#define PARAM_T 2
#define PARAM_P 3

#define FOR_EACH_PARAM(I) for(I = 0; I < 4; I++)

#define a(X) value[PARAM_A][X]
#define b(X) value[PARAM_B][X]
#define t(X) value[PARAM_T][X]
#define p(X) value[PARAM_P][X]

/* store which iteration each param is being calculated */
int iter[4] = {1, 1, 1, 1};

/* store all parameters values */
mpf_t *value[4];

pthread_t thread[4];
/* mutexes to protect simultaneous access to 'iter' variable */
pthread_mutex_t mutex[4];

/* routines to calculate the parameters a, b, t and p */
void *calc_param_a(void *);
void *calc_param_b(void *);
void *calc_param_t(void *);
void *calc_param_p(void *);

void * (*func[4])(void*) =
    {calc_param_a, calc_param_b, calc_param_t, calc_param_p};

int main(int argc, char *argv[])
{
    int n_iter = DEFAULT_N_ITERATIONS;
    size_t prec = DEFAULT_PRECISION;

    /* First argument overrides the default number of iterations */
    if(argc > 1)
        sscanf(argv[1], "%d", &n_iter);

    /* Second argument overrides the default precision in bits */
    if(argc > 2)
        sscanf(argv[2], "%zu", &prec);

    mpf_set_default_prec(prec);

    int i;
    FOR_EACH_PARAM(i)
        value[i] = (mpf_t *) malloc(n_iter * sizeof(mpf_t));

    /* Initialize all variables in iteration 0 */
    mpf_init_set_d(a(0), 1.0);

    mpf_init_set_ui(b(0), 2);
    mpf_sqrt(b(0), b(0));
    mpf_div(b(0), a(0), b(0));

    mpf_init_set_d(t(0), 0.25);

    mpf_init_set_d(p(0), 1.0);

    /* Initialize mutexes */
    FOR_EACH_PARAM(i)
        pthread_mutex_init(&mutex[i], NULL);

    /* Initialize threads */
    FOR_EACH_PARAM(i)
        pthread_create(&thread[i], NULL, func[i], &n_iter);

    /* Calculate pi \o/ */
    mpf_t *a, *b, *t;
    mpf_t pi;

    mpf_init(pi);

    /* Wait for parameters a and b */
    pthread_join(thread[PARAM_A], (void **) &a);
    pthread_join(thread[PARAM_B], (void **) &b);

    /* Operate them while t is not ready */
    mpf_add(pi, *a, *b);
    mpf_mul(pi, pi, pi);

    /* Wait for parameter t */
    pthread_join(thread[PARAM_T], (void **) &t);

    /* Finish to calculate pi */
    mpf_mul_ui(*t, *t, 4);
    mpf_div(pi, pi, *t);

    fprintf(stderr, "pi\n");

    mp_exp_t exponent;
    char *result = mpf_get_str(NULL, &exponent, 10, 0, pi);

    printf("%s\n", result + exponent);

    free(result);

    /* Clear all the mess :D */
    int j;
    FOR_EACH_PARAM(i)
    {
        for(j = 0; j < n_iter; j++)
            mpf_clear(value[i][j]);

        free(value[i]);
    }

    return 0;
}

void *calc_param_a(void *arg)
{
    int n_iter;
    /* We know that b(0) was already calculated */
    int b_iter = 1;

    n_iter = * (int *) arg;

    int i;
    while(iter[PARAM_A] < n_iter)
    {
        i = iter[PARAM_A];

        /* Verify if b(i-1) is ready */
        while(b_iter < i)
        {
            pthread_mutex_lock(&mutex[PARAM_B]);
            b_iter = iter[PARAM_B];
            pthread_mutex_unlock(&mutex[PARAM_B]);
        }

        mpf_init(a(i));

        mpf_add(a(i), a(i-1), b(i-1));
        mpf_div_ui(a(i), a(i), 2);

        /* Lock all threads that depends on a */
        pthread_mutex_lock(&mutex[PARAM_A]);

        /* Update iteration */
        iter[PARAM_A]++;

        /* Unlock threads that depends on a */
        pthread_mutex_unlock(&mutex[PARAM_A]);

        fprintf(stderr, "a(%d)\n", i);
    }

    return &a(i);
}

void *calc_param_b(void *arg)
{
    int n_iter;
    /* We know that a(0) was already calculated */
    int a_iter = 1;

    n_iter = * (int *) arg;

    int i;
    while(iter[PARAM_B] < n_iter)
    {
        i = iter[PARAM_B];

        /* Verify if a(i-1) is ready */
        while(a_iter < i)
        {
            pthread_mutex_lock(&mutex[PARAM_A]);
            a_iter = iter[PARAM_A];
            pthread_mutex_unlock(&mutex[PARAM_A]);
        }

        mpf_init(b(i));

        mpf_mul(b(i), a(i-1), b(i-1));
        mpf_sqrt(b(i), b(i));

        /* Lock all threads that depends on b */
        pthread_mutex_lock(&mutex[PARAM_B]);

        /* Update iteration */
        iter[PARAM_B]++;

        /* Unlock threads that depends on b */
        pthread_mutex_unlock(&mutex[PARAM_B]);

        fprintf(stderr, "b(%d)\n", i);
    }

    return &b(i);
}

void *calc_param_t(void *arg)
{
    int n_iter;
    /* We know that a(0) and p(0) were already calculated */
    int a_iter = 1;
    int p_iter = 1;

    n_iter = * (int *) arg;

    int i;
    while(iter[PARAM_T] < n_iter)
    {
        i = iter[PARAM_T];

        /* Verify if a(i) is ready */
        while(a_iter < i + 1)
        {
            pthread_mutex_lock(&mutex[PARAM_A]);
            a_iter = iter[PARAM_A];
            pthread_mutex_unlock(&mutex[PARAM_A]);
        }

        mpf_init(t(i));

        /* Partially calculate t(i) */
        mpf_sub(t(i), a(i-1), a(i));
        mpf_mul(t(i), t(i), t(i));

        /* Verify if p(i-1) is ready */
        while(p_iter < i)
        {
            pthread_mutex_lock(&mutex[PARAM_P]);
            p_iter = iter[PARAM_P];
            pthread_mutex_unlock(&mutex[PARAM_P]);
        }

        mpf_mul(t(i), t(i), p(i-1));
        mpf_sub(t(i), t(i-1), t(i));

        /* Lock all threads that depends on t */
        pthread_mutex_lock(&mutex[PARAM_T]);

        /* Update iteration */
        iter[PARAM_T]++;

        /* Unlock threads that depends on t */
        pthread_mutex_unlock(&mutex[PARAM_T]);

        fprintf(stderr, "t(%d)\n", i);
    }

    return &t(i);
}

void *calc_param_p(void *arg)
{
    int n_iter;

    n_iter = * (int *) arg;

    int i;
    while(iter[PARAM_P] < n_iter)
    {
        i = iter[PARAM_P];

        mpf_init(p(i));

        mpf_mul_ui(p(i), p(i-1), 2);

        /* Lock all threads that depends on p */
        pthread_mutex_lock(&mutex[PARAM_P]);

        /* Update iteration */
        iter[PARAM_P]++;

        /* Unlock all threads that depends on t */
        pthread_mutex_unlock(&mutex[PARAM_P]);

        fprintf(stderr, "p(%d)\n", i);
    }

    /* It's not necessary return p because we don't use it */
}
